Code
if(!require(sf)) install.packages("sf")Loading required package: sf
Linking to GEOS 3.11.2, GDAL 3.6.2, PROJ 9.2.0; sf_use_s2() is TRUE
So far, we have largely focused on the non-geographic data set, mtcars. As someone engaged with a geographic discipline, it is highly likely that you will be engaging geographical data frames. Such data sets have the advantage of being able to produce cartographic representations - maps!
Principally, R is a statistical script-based programming language without mapping capabilities. To produce maps in R you’ll need to read in a package called sf. sf is an abbreviation of Simple Features (Pebesma and Bivand 2023) that handles a data-frame objects with simplified code-lines including geographic formats.
Maps are initially complex to make in R but offer another way of describing data with geographic properties.
Geographic data come in variety of formats but essentially categories as either being vector-based or rasterised. For now, we’ll focus on manipulating a vector data frame containing polygons. First, to load sf if you’ll need to install package and then load its libraries into your R session. A way doing this with minimal as is as follows:
Loading required package: sf
Linking to GEOS 3.11.2, GDAL 3.6.2, PROJ 9.2.0; sf_use_s2() is TRUE
When running a line of code similar that above, in a script, you may need to run it twice for the package to install and then to load. To avoid this, you can also use the following approaches:
If you have already have the package installed you may want to comment-out the first line by placing a hashtag before it:
An example of geographic data can Tenure types recorded in Bristol, at the Lower Super Output Area (LSOA) level: https://raw.githubusercontent.com/Richtea84/I2Q-files/main/Bristol_Tenure.csv
It is possible to load the CSV into R as an object using the weblink above in using the read.csv(...) function. Let’s call the Object ‘Bristol_data’.
You can preview the results using the head(...) function. For example:
X2021.super.output.area...lower.layer Total..All.households Owned
1 E01014601 : Bristol 001A 746 474
2 E01014602 : Bristol 001B 827 353
3 E01014603 : Bristol 001C 1068 604
4 E01014605 : Bristol 001E 835 237
5 E01032516 : Bristol 001G 859 578
Shared.ownership Social.rented Private.rented Lives.rent.free
1 5 187 77 3
2 6 373 95 0
3 27 195 242 0
4 1 473 120 4
5 7 139 132 3
Owns.with.a.mortgage.or.loan.or.shared.ownership
1 247
2 196
3 371
4 137
5 252
Private.rented.or.lives.rent.free
1 80
2 95
3 242
4 124
5 135
As discussed in previous sections, there are ways in which you can change the column names - for now, we’ll leave them as they are. The data we have, is a ‘data frame’.
As we are looking at 2021 LSOAs for Bristol, we can access the corresponding geographic boundaries from UK data service: https://borders.ukdataservice.ac.uk/
Once there do the following:

You may encounter further prompts related to the selection of the data download format - ensure you pick the ‘Download features in Shapefile format as Zip file’ for your selected data set when prompted.
At this point, it is important that a relevant working directory has been set up. The directory will contain both you CSV and the zipped data that you’ve downloaded. Next, boundary data is unzipped and subsequently loaded into R as a new object (called Bristol_LSOAs) using the read_sf(...) function from the sf package:
For the shapefile (.shp) to load correctly, ensure that all of the other dependencies associated with the shapefile are present. Here, these are:
Under some circumstances, other dependencies may be present such as:
Whatever comes across in the download, ensure that it is present in the folder from which you are reading the .shp file.
The process of joining CSV data to the mappable data uses the merge(...) function in exactly the same way as before. An important detail is that the mappable data is always the first object to be called in.
First, let’s change the name of the field containing the containing the LSOA code reference (beginning ‘E01’).
The ‘code’ field we’ve just created contains too much information, making it impossible to match to Bristol_LSOAs data. The LSOA code is found in the first 9 characters. To shorten the contents here, we can use the substr(...) function:
x = is the same fieldstart = specifies the start of the crop - here, it is the first characterstop = specifies the end of the crop - here, it is the ninth characterWith the fields properly prepared, we are able to conduct the join using the merge(...) function. Again, note that the geographic data comes first (x=).
With sf loaded, it is possible to preview the result of the join using the plot function. Let’s examine the first 6 plots using a slicing:
Boundary of the Bristol featured in the LSOA data features the Avonmouth Extension (Figure 1) that protudes out into the sea due the region’s in relation to Bristol’s maritime history that saw shipping trade attributed to the large tidal range of the River Avon that enabled Bristol’s Floating Harbour; shipping traffic falls under the jurisdiction of Bristol’s port authority within the extension. However, we are only interested with the land administration component, where we will need to clip the region accordingly.
The first step is to ‘clip’ the extension off so that we are only looking at the land. This is achieved by downloading the boundaries for Bristol’s wards. The Shapefile version of the data should be downloaded and loaded following a similar process to that described in Section 0.2.2.
Let’s load the shapefile as an object called ‘Bristol_wards’ (ensuring all of the dependencies are present):
We can plot one of the variables using the plot(...) function:
This good, but let’s merge all of the separate polygons into a single landmass. The sf package features a useful function that dissolves boundaries called st_union(...):
Now that we have the dissolved (or unionised) layer, we are in a position to clip off the Avonmouth Extension. Here’s we’ll use the st_clip(...) function and use the Bristol_Mass object as the mask.
x = The data that needs to be cropped (Bristol_merged)y = The cropping backgroundWarning: attribute variables are assumed to be spatially constant throughout
all geometries

At this point, we are examining continuous data where the blue tones indicate lower observation frequencies, while the yellow tones are the higher frequencies.
Inspecting social rentals, we can detect areas where social rentals (i.e. council tenancies) are highest - in the middle and towards the bottom. However, it is difficult to ascertain scale and location, where this graphic lacks cartographic elements and cannot be considered a ‘map’ at this stage. To introduce cartographic elements another package needs to be introduced, tmap (Tennekes 2018). Loading this package follows a similar process to that described above:
The parsing for tmap differs from the default R coding structure and closely aligned to an extended graphical suite called ggplot that features in other programming languages such as Python. Here’s the code for producing a basic preview for our clipped data with tmap is as follows:

tmaptm_shape(...) introduces the geographic (mappable data)tm_polygons(...) specifies the geometry and its aesthetic properties.These are linked together with + symbols.
To make this graphic cartographic, we need to include the following elements:
tm_compass(...)tm_scale_bar(...)tm_graticules(...)Adding these on can be achieved by using further + symbols in your code. An example is shown in Figure 5.

As with the plot(...) function, colours and elements can be manipulated. You can learn more about how to customise the look and feel of your map by inserting ? in front of the graphical element that you want to manipulate.
An example of how arguments can be applied to improve the map above is as follows:
tm_shape(shp = Bristol_merge_2)+
tm_polygons(col = "Social.rented", title = "Social rented", border.alpha = 0)+
tm_compass(position = c("right", "top"), color.dark = "navy",type = "4star", text.color = "navy")+
tm_scale_bar(position = c("left", "bottom"), color.dark = "navy",text.color = "navy", breaks = c(0,2,4,6))+
tm_graticules(col = "navy")+
tm_layout(legend.outside = TRUE, legend.position = c("right","center"),
bg.color = "grey", frame = FALSE)
In addition to the placement, the map’s default colour scheme can also be manipulated. Using either a default Brewer palette, importing a palette library, or by creating your own customised palette. Here’s an example of a map produced using a custom colour palette (Figure 7)
custom_pal <- c("red","grey20")
tm_shape(shp = Bristol_merge_2)+
tm_polygons(col = "Social.rented", title = "Social rented", border.alpha = 0, palette = rev(custom_pal), n = 12)+
tm_compass(position = c("right", "top"), color.dark = "navy",type = "4star", text.color = "white")+
tm_scale_bar(position = c("left", "bottom"), color.dark = "navy",text.color = "white", breaks = c(0,2,4,6))+
tm_graticules(col = "navy")+
tm_layout(legend.outside = TRUE, legend.position = c("right","center"),
bg.color = "grey", frame = FALSE)
A critical step in the production of this graphic is the use of palette = argument that points to the custom colour palette. Within this argument the rev(...) function is applied to reverse the colour palette.
While someone with a good understanding of geometry can figure out the precise location of the featured regions, for most the inclusion of a base map is far more practical. This can be added using a package called rosm (Dunnington 2022) that has access to Open Street Map map tiles stored in a raster format. We’ll discuss rasters in another section.
As Open Street Map is projected to the Mercator projection (same as Google Maps), we need to change the projection of our map data (Bristol_merge_2) so that it conforms. The EPSG code we need is 4326 This can be achieved using the st_transform(...) from the sf package. Here is how it is applied in the creation of a new object simply called Bristol_reproj:
With the application of the Mercator projection assigned to our data, we can now retrieve a relevant base map. For this, we need the osm.raster(...) function from rosm and the st_bbox(...) function from sf to define the bounding box (or rectangular extents) of our data.
To add the base map, a raster image, we use tm_shape(...) in combination with the tm_rgb(...) function from tmap. Here’s our base map:
As with other elements in tmap the basemap can be manipulated. Let’s make it monochrome by lowering the saturation = argument to 0:

Adding our mapped data onto of this is this base map is achieved by adding the map code immediately after the base map. Transparency on for out polygon layer is achieved by introducing the alpha = argument, the selected value is 0.45:
tm_shape(my_basemap)+
tm_rgb(saturation = 0)+
tm_shape(shp = Bristol_merge_2)+
tm_polygons(col = "Social.rented", title = "Social rented", border.alpha = 0, palette = rev(custom_pal), n = 12, alpha = 0.45)+
tm_compass(position = c("right", "top"), color.dark = "navy",type = "4star", text.color = "white")+
tm_scale_bar(position = c("left", "bottom"), color.dark = "navy",text.color = "white", breaks = c(0,2,4,6))+
tm_graticules(col = "navy")+
tm_layout(legend.outside = TRUE, legend.position = c("right","center"),
bg.color = "grey", frame = FALSE)
As seen in Figure 10, we are able to contextualise data we’ve imported by seeing the names of nearby areas. Looking at our web browswers it is possible to identify that Hartcliffe and areas close to Easton, Redfield, and St Philips are have high occurrences of social rented housing.
A helpful feature in tmap is its ability to produce interactive web maps. These are essentially rendered using a javascript package known as leaftlet. Creating an interactive map in tmap is achieved by using the following command before your main code tmap_mode("view"). Here’s the result when applied to our map above:
tmap mode set to interactive viewing
#tm_shape(my_basemap)+
#tm_rgb(saturation = 0)+
tm_shape(shp = Bristol_merge_2)+
tm_polygons(col = "Social.rented", title = "Social rented", border.alpha = 0, palette = rev(custom_pal), n = 12, alpha = 0.45)+
tm_compass(position = c("right", "top"), color.dark = "navy",type = "4star", text.color = "white")+
tm_scale_bar(position = c("left", "bottom"), color.dark = "navy",text.color = "white", breaks = c(0,2,4,6))+
tm_graticules(col = "navy")+
tm_layout(legend.outside = TRUE, legend.position = c("right","center"),
bg.color = "grey", frame = FALSE)Compass not supported in view mode.
legend.postion is used for plot mode. Use view.legend.position in tm_view to set the legend position in view mode.
Warning: In view mode, scale bar breaks are ignored.
A crucial concern with this mode of mapping is the loss of cartographic elements including the graticules, and compass - this results in a number of warnings when our raw code is introduced. Since several base maps are supplied by default, we no longer have to generate a base map. To change the mode back, we use the following line of code: